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ABSTRACT 



Inversions of rotational splittings have shown that 
the surface layers and the so-called solar tachocline 
at the base of the convection zone are regions in which 
high radial gradients of the rotation rate occur. The 
usual regularization methods tend to smooth out ev- 
ery high gradients in the solution and may not be ap- 
propriate for the study of a zone like the tachocline. 
In this paper we use nonlinear regularization methods 
that are developed for edge-preserving regularization 
in computed imaging (e.g. Blanc-Feraud et al. 1995) 
and we apply them in the helioscismic context of ro- 
tational inversions. 



1. INTRODUCTION 



Moscgaard 1996] ). An investigation of the possibil- 
lty to use the elaborate nonlinear techniques, devel- 
oped for edge-preserving regularization in computed 
imaging, in hclioseismic context is being developed by 



Corbard et al. (1998b). Here we present preliminary 



results obtained by this method for the tachocline. 

Section |] briefly recalls the relation between the so- 
lar internal rotation and the helioseismic measured 
frequency splittings and presents the corresponding 
discretized inverse problem. We discuss in Section p] 
the non linear approach of regularization in inverse 
techniques and the computational aspects. In Sec- 
tion the choice of the regularizing parameters for 
the particular case of the solar rotation inversion and 
the preliminary results obtained with LOWL (Tom- 
czyck et al. 1995, 2 years of data) and MDI (144 
days, Schou et al. 1998) data are presented. 



The existence of high gradients in the solar rotation 
profile near the surface and at the base of t he convec- 



tion zone i n the so-called solar tachocline ( Spiegel & 
Zahn 199*2] ) has been reveale d by the inversions ot the 



rota tionaTsplittings (see e.g. Thompson et al. (1996) 
and Schou et al. (1998) for the last results J. The 
tachocline represents a thin zone where the differ- 
ential rotation of the convection zone becomes rigid 
in the radiative interior. It is thought to be the 
place from where the solar dynamo originates and 
its precise structure is an important test for angu- 
lar momentum transport theories. More precisely, 
the thickness of the tachocline can be related to the 
horizontal component of the turbulent viscosity and 
may be used as an important observa tional constraint 
on the properties of the turbulence ( Spiegel fc Zahn 



1992 



Gough fc Sekii 1998; Elliot 1997 J 



Several works have already been per formed to in 



fcr the fine structure of the ta chocli ne fCharbonncau 

Bgngg |Cor 



et al. 



1998j: jKosoyichev 1996|. 1998: 
bard et al iy98a| ; |Antia~cT~al. 1998[ ' using both tor- 
ward analysis and inverse techniques. For the in- 
verse approach, it may be interesting to change the 
global constraint which tends to smooth out every 
high gradients in the solution and to find a way to 
preserve such zones in the inversion process. A first 



attempt in this di rection has been carried out (Cor- 



bard et al. 1998a) by using a nonlinear regulariza- 



lon term through the PP-TSVD method (Hansen & 



2. FORMULATION OF THE PROBLEM 



In this paper we consider the ID problem of inferring 



the solar equatorial rotation profile Q e 



fil 



from sectoral splittings (i.e. 
1984; [Antia et al. 1996j) 



7-^90° 



I) (e.g. Duvall et al 



Vnll — VnlO 

I 



K ni {r) n eq (r) dr, (1) 



where /, n, m are respectively the degree, the ra- 
dial order and the azimuthal order, r is the solar 
radius and K n i (r) are the so-called rotational kernels 
that have been calculated for each mode from a solar 



model taken from Morel et al. (1997) 



This approximation of the 2D integral equation which 
relates the internal rotation to th e splittings, is valid 



only for high degree modes (e.g. Corbard 1997) but 
the influence of low degree modes on the determina- 
tion of the tachocline and upper layers is thought to 
be small. Moreover the rotation is known to be rigid 
or quasi rigid in the radiative interior. 

We search a solution H(r) as a piecewise linear func- 



2 



tion of the radius by setting: 

^( r ) = y^^p^pi 7 ") ^ = (^p)p=i,A r P (2) 
P =i 

where Vv( r )j P — l>-^p are piecewise straight lines 
(iVp = 100 in this work) between fixed break points 
distributed a ccording to the dens ity of turning points 
of modes (cf. JCorbard et al. 1997 ) . The discretization 
of Equation R\ leads to the matrix equation: 



well constrained by the data. In this work we con- 
sider only smoothing term with first derivative. 

For a general (^-function, one can write the criterion 
Equation |in a discretized form: 



J(17)= X 2 (17) + A 2 J 2 (ft), 



(8) 



where J2(f2) represents the discretized regularization 
term defined by: 



w = m, 



(3) 



where W = (Wi/<Ti)i=i t N is the vector of the N ob- 
served frequency splittings Wi weighted by the stan- 
dard deviation Oi given by observers for each mode 
i = (n, I). No correlation between the different modes 
is assumed. The matrix R is defined by: 



R={R 



Rip — 



1 



R-. 



Ki(r)ip p (r)dr (4) 



An inverse method should lead to a solution that is 
able to produce a good fit of the data. We define the 
goodness of the fit in chi-square sense by the x 2 value 
obtained for any solution f2(r): 



■•P 



dr 



dr 



N p -l 

£ 

P =i 



Cpip 



(I 



m\ 



. ( 9 ) 

In this equation (c p ) p= i.n p -i represent the weights 
used for the integration rule, L is a discrete approx- 
imation of the first derivative operator, and |if2| is 

the absolute value of the p-component of t he vector 
Lfl. The expre ssion of c p and L are given in Corbard 
et al. (1998b) for the simple case of the polynomial 



expansion Equation || 

The minimization of the criterion J(Q) leads to the 
following Euler equations (discretized form): 



! (^)) = E 



nl 



Wi- f* Ki(r)n{r)dr 



, (5) 



VJ(n) = (r t r + x 2 L T B(n)L)n = r t w 

( 10 ) 

B is a diagonal matrix. Its elements depend on the 
gradient of the solution at each grid point: 



which can be written in the discretized form: 



x 2 (n) = \\Rn-w\\l 



(6) 



B = diag(b p ) with b p = c p x 



<P (\L%) 

J (11) 



2\m\ 



NON LINEAR REGULARIZATION 



3.1. Euler equations 

The inverse integral problem is an ill-posed problem 
and the minimization of only the x 2 value generally 
leads to oscillatory solutions that are not 'physically 
acceptable'. So, a regularization technique must be 
used in the minimization process. A large class of 
these techniques can be expressed in the general form 
of the minimization of a criterion J over the unknown 
solution Vl (r): 



J(n(r))= X 2 (n(r)) + \< 



n., : 



d q n{r) 



dri 



dr, (7) 



The so-called trade-off parameter A is chosen so that 
it establishes a balance between the goodness of the 
fit of the data and the constraint introduced on the 
solution (cf . Section [I]) . The order q of the deriva- 
tive is usually taken equal to one or two. The two 
choices can lead to similar results with the appro- 
priate choice A in the domains where the solution is 



Two choices for the (/^-function lead to well known 
regularization strategies: 



ip(t) = t 2 corresponds to the traditional 
Tikhonov approach with first derivative whereas 

tp(t) = t is known as the Total Variation (TV) reg- 
ularization method (e.g. Acar fc Vogel 1994 ). It 
has been shown that this regularization method 
is able to recover piecewise smooth sol utions with 
steep gradients ( Vogel fc Oman 1996| ) . 



The use of a general nonquadratic ^-function will 
lead to a nonlinear problem which requires an appro- 
priate iterative method to be solved. 



3.2. Properties of the weight function 



2t 



we 



From the Euler equation (Equations |10| and 

can see that the function <p (t)/2t acts as a weight 
function in the smoothing process: at each grid point 
the gradient of the solution is used as an argument of 



3 



this function in order to set locally more or less reg- 
ularization. This suggests an iterative process where 
th e gradient of the solution at a given step is used 



for the computation of the regularization term at the 
next step. Three properties of the weighting function 

if (t)/2t are required to obtain a s atisfactory solution 
and to preserve high gradients (Charbonnier et al 
19Q: 



1. no smoothing for high gradients: 

lim^=0 

2. Tikhonov smoothing for low gradients: 



(12) 



< lim = M < oo (13) 



to Equation ll| with addition of Gaussian noise with 
a standard deviation taken for each mode from the 



formal error given in observational data (cf. Corbard 
et al. 1998a| ) . A second set of artificial data with the 



same rotational law and standard deviations divided 
by a/10 has also been used. 

At each stepof ARTUR algorithm the linear system 
(Equation [lj) has been solved using an iterative con- 
jugate gradient method using f2 fc_1 as starting point. 
This leads to a very fast algorithm where the num- 
ber of conjugate gradient iterations needed to solve 
the linear system decreases at each ARTUR step (i.e. 
as k increases). The algorithm is stopped when the 
norm of the relative difference between two solutions 
at two successive steps is below 10 -6 i.e.: 



\n k - n 



k-li 



\n k \ 



< 10 



-6 



(16) 



3. 



<P (*) 
2t 



strictly decreasing to avoid instabilities. 

(14) 



Either c onvex or non-convex ^-func tion s imav be cho- 
sen (se e |Charbonnier et al. (1997) and Tcboul et al, 
(1998)| for examples in both cases). A non-convex 
function may be more suited for the search of high 
gradients. But this choice leads to some numerical 
difficulties and instabilities related to the existence 
of local minima and may induce a high sensibility 
of the inverse process to the regularization parame- 
ters. At the opposite the choice of a convex function 
avoids these numerical problems and is more suitable 
for relatively smooth transition (Blanc-Feraud et al 
19951). 



3.3. The iterative algorithm: ARTUR 



Following Charbonnier et al. (1997) the inversion us- 
ing non linear regularizing criterion can be solved by 
an iterative scheme named ARTUR (Algebraic Re- 
construction Technic Using Regularization) that is 
easy to implement: at each step k we calculate the 
regularization term using the derivative of the previ- 
ous estimate £l k ~ 1 and we simply compute the new 
estimate fl k by solving the linear system: 



R T R + X 2 L T B(n k - 1 )L)Q k = R T W. (15) 



For a convex ^-function, the convergen ce of this so- 
called algor ithm has been established ( |Charbonnier 
et al. 1997|). This is therefore an adaptative regu- 



larization method which uses the information on the 
derivative of the solution obtained at each step in or- 
der to improve the regularization at the next step. 
This requires an initial guess fl° for the solution but 
we will show in the next section that a constant so- 
lution can ever be used as the starting point. 

An example of artificial discontinuous rotation has 
been used to test the algorithm. The corresponding 
rotational splittings have been computed according 



The results are given in Figures @ and | which show 
examples of ARTUR steps (upper window) . It is seen 
how the ARTUR solution, starting from a smoothed 
Tikhonov solution, becomes steeper at each step. 
Comparison between Figures [0 and @ shows the effect 
of the errors in the data which leads to a smoothing 
of the edges of the discontinuity. 





ARTUR Steps 
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Figure 1. Solutions obtained by inverting splittings com- 
puted from a discontinuous one dimension rotation pro- 
file (full line) for the same mode set as in LOWL data 
and by adding some 'realistic' noise (see text). The stan- 
dard Tikhonov solution is given for two different auto- 
matic choices of the regularizing parameter. The succes- 
sive steps of ARTUR algorithm are given in the upper left 
window whereas the final step is shown on the main plot. 
The choice of the regularizing function and parameters 
for ARTUR algorithm are those discussed in Section u\ 
The solutions are plotted without error bars for clarity. 



APPLICATION TO THE SOLAR ROTATION 
INVERSION 



In the particular case of the solar tachocline, the un- 
certainty on the widt h of the transition zone i s still 
large (see Table 2 of Corbard et al. (1998a) for a 



4 




Initia rotation rate 
Tikhonov + L— curve 
Tikhonov + GCV 
ARTUR (final step) 



0.7 
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Figure 2. The same as Figure ^ but computed for a lower 
level of noise (standard deviations divided by y/lO). Com- 
parison between Figures ana! @ shows the smoothing ef- 
fect of the data noise level for the three methods. 



summary of some previous works). Here we apply 
the non linear regularizing method to the LOWL and 
MDI sectoral splittings. 



4.1. Choice of the (/^-function 



4.2. The choice of A 

If the initial guess f2° is a constant function, then, 
according to Equations Ol and |l9L M — 1 and the 
solution at the first ARTUR step corresponds to a 
Tikhonov solution wi th A as regularizing pa rameter. 
It has been shown in |Corbard et al. (1998a) that the 
Generalized Cross Validation (GCV) strategy leads 
systematicall y to a less sm oothed solution than the 
L-curve one ( Hansen 1992[ ) (Xhcurve — 100 * Xgcv in 
that work) and therefore is more suited to the study 
of the tachocline. Nevertheless, this choice leads to 
spurious oscillations below and above the tachocline 
(see Figures |, 0). As ARTUR algorithm will tend to 
enhance the nigh gradients found at the first step it 
is important to start with a solution smooth enough 
to avoid spurious oscillations with high gradients. At 
the opposite, the L-curve choice leads often to a so- 
lution which is too smooth and does not allow to 
exhibit the expected high gradients during the it- 
erations. As the optimal choice of this parameter 
strongly depends on the level of noise included in the 
data (cf. Figures 1 and 2), it is important to define 
an automatic choice of this parameter so that we use 
the same strategy for different datasets or for differ- 
ent realizations of the noise. It consists in taking an 
intermediate value between the values of \hcurve and 
Xgcv f° r Tikhonov inversion. This has been used 
here to fix A in ARTUR algorithm. 



According to the previous discussion, we have cho- 
sen to consider a convex regularizing ^-function as 



4.3. The choice of 5 



Charbonnier et al. (1997) 



ip(t) = 2\Jt 2 + 1-2 



(17) 



This function is close to the absolute value function 
used in TV regularization but has a quadratic be- 
havior near in order to smooth zones with small 
gradient moduli in addition to the linear behavior 
near infinity that preserves high gradient zones. 

In addition to the usual regularizing parameter A al- 
ready introduced, we have included in the ^-function 
a parameter which scales the derivative of the solu- 
tion. It allows to adapt the domain of gradient mod- 
ulus where the solution is very weakly smoothed to 
the specific problem we consider. 





150 200 
t=|dQ/dr| 



J(fi(r))= X 2 (^(r))+A 2 



dfl(r) 



dr 



dr 
( 18 ) 

This simply leads to replace A by A = A/5 in Euler 
Equation Hfl and to use 



Figure 3. The full line shows the first derivative of a first 
step solution (cf. Figure J^) in ARTUR algorithm as a 
function of the fractional solar radius. For each value of 
the gradient, the dashed line gives the weight that will be 
given locally to the regularizing term at the second step of 
ARTUR algorithm. The dotted line indicates the gradient 
upon which the local regularization will be more than 50% 
smaller at the second step compared to the first step. 



tp(t)/2t=l/<Jl + (t/S) 



(19) 



where t = \LQ\ as weighting function in Equa- 
tion |n]. 

Therefore we have to define a strategy to choose the 
two parameters A and S. 



The parameter S is introduced to adapt the shape of 
the weighting function to the gradient that we search 
to detect . Its value is chosen by looking at the deriva- 
tive of the solution at the first iteration step. We have 
chosen for simplicity to keep this parameter constant 
during the iterations. Figure ^ shows as an example 
(full line) the first derivative of solution obtained at 
the first step by inverting artificial splittings which 



5 



1 


IT 
if 
ll 


,r \ 




Jl 
A J 





0.8 0.9 1 







/N - 


- 


it 

V 




ft 

i 1 

y i 
s 





0.5 0.6 



0.7 



Figure 4- Solar equatorial rotation obtained by invert- 
ing LOWL data. Tikhonov solution is shown by the full 
line with error bars. The dashed line represents the final 
ARTUR step. 



have been computed for the discontinuous rotation 
law presented in the previous section (cf. Figure |l|). 
The largest peak corresponds to the rapid variation 
in the tachocline, the smaller ones to spurious oscilla- 
tions in the solution. The weighting function (Equa- 
tion |l^) is shown in dashed line for 8 = 100. Ac- 
cording to Figure the choice of 8 = 100 leads to 
regularize 50% less at the second step in that zones 
where the gradient of the first step solution is above 
- 175 nHz/i? Q . 

According to the previous results on solar rotation in- 
versions, the width of the tachocline does not exceed 
0.1 solar radius. This represents also the resolution 
reached at the tachocline localization (~ 0.69i?o) 
with a Tikhonov method using a regularizing param- 
eter chosen near the corner of the L-curve. Further- 
more we have a good estimate of the difference be- 
tween the rotatio n rate above and below the transi- 
tion (~ 30nHz in |Corbard et al. (1998a) ) . Therefore 
we can estimate an order of magnitude of 3OOnHz/i?0 
for the maximum gradient obtained at the first iter- 
ation of ARTUR process. This corresponds to the 
value in Figure 0. At the second step we want to 
preserve only high gradients i.e. to regularize less in 
that zones where high gradients have already been 
found at the first step. With our a priori knowledge 
of the maximunigradient at the first step (~ 300nHz, 
see also Figure ||), the choice of 8 = 100 sounds rea- 
sonable in the sense that it will tend to decrease the 
regularization especially in the tachocline. A smaller 
value would enhance the secondary peaks that may 
be induced by the data noise. 



4.3.1. Results for the tachocline 



Figures [l] and ^ show that ARTUR algorithm leads 
always to better results than Tikhonov inversion 
without 'local deconvolution' in the case of a discon- 
tinuous rotation profile . Furthermore, it is shown in 
Corbard et al. (1998b)| that it can also lead to good 
results tor simulated tachocline widths between 0.02 
and O.O8i? . 



Figure 5. Same as Figure 5 for MDI data 



Preliminary results have been obtained by inverting 
the sectoral splittings given by the 2 years LOWL 
data and the 144 days MDI data. The results of 
Tikhonov inversion and the last step of ARTUR al- 
gorithm are plotted in Figures Q and pi Tikhonov so- 
lutions are shown with error bars at each grid points 
which are deduced from the propagation of data noise 
through this linear process. ARTUR algorithm is 
nonlinear and therefore we can not compute formal 
errors in the same way. 

In order to have an estimate of the uncertainties on 
the tachocline parameters derived by ARTUR algo- 
rithm, Monte-Carlo simulations with artificial split- 
tings have been performed by Corbard et al. (1998b). 
The rotation profiles were taken as erf functions with 
different 'initial widths'. Each 'inferred width' is the 
mean value of the results obtained by fitting directly 
the solutions by an erf function, for 500 realizations 
of input errors. Error bars are estimated from a 
68.3% confidence interval. From this study it seems 
that an uncertainty of ±0.02i? Q on the tachocline 
width can be reached with both ARTUR algorithm 
and Tikhonov inversion with a 'local deconvolution' 
using the averaging kernels computed at the center 
of the transition (see Corbard et al. 1998b). 

The width of the tachocline estimated from LOWL 
data and by using ARTUR algorithm is 0.05±0.02i? Q 
in good agreement with t he value of 0.05 ± 0.03i ?^ 
found for the same data in Corbard et al. (1998a)| by 
studying systematically the ettcct ot regularization on 
the determination of tachocline parameters for three 
inverse methods. A preliminary study of MDI data 
leads to a slightly larger width of 0.08i?o but Monte- 
Carlo simulation have not yet been performed that 
allows to give an estimate of the uncertainty on this 
value. For both dataset the widths obtained from 
Tikhonov inversions after 'local deconvolution' using 
averaging kernels are always larger (~ O.li?©) than 
the estimates obtained with ARTUR algorithm. 



We note that these estimates of the tachocline width 
have been obtained by fitting the solutions by an erf 
function between 0.4 and 0.8i?o and this may be bet- 
ter adapted to the shape of LOWL solution (Figure 0) 
which presents a step between 0.75 and O.8i?0. This 
step is not found with MDI data and this may ex- 
plain the larger width found by our fit with these 
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data. One possibility for future work is to change 
our fitting function i.e. to change our definition of 
the tachocline width. However we have also to ex- 
plain the different behavior of the two solutions near 
O.8i? and to look if this remains with longer time 
series of MDI experiment. 



ACKNOWLEDGMENTS 



We acknowledge SOI team and S. Tomczyk for al- 
lowing the use of MDI and LOWL data. SoHO is 
a project of international cooperation between ESA 
and NASA. This work has been performed using the 
computing facilities provided by the program "Sim- 
ulations Interactives et Visualisation en Astronomie 
et Mecanique" (SIVAM, OCA, Nice) and by the "In- 
stitut du Developpement et des Ressources en Infor- 
matique Scientifique" (IDRIS, Orsay). T. Corbard 
thanks the conference organizers for financial sup- 
port. 



REFERENCES 

Acar, R., Vogel, OR. 1994, Inverse Problems 10(6), 
1217 

Antia, H.M., Basu, S., Chitre, S.M. 1998, MNRAS, 
submitted 

Antia, H.M., Chitre, S.M., Thompson, M.J. 1996, 
A&A 308, 656 

Basu, S. 1997, MNRAS 288, 572 

Blanc-Feraud, L., Charbonnier, P., Aubert, G., Bar- 
laud, M. 1995, In: IEEE Proceedings of the 
2nd International Conference of Image Processing, 
Washington DC, USA, p. 474 

Charbonneau, P., Christensen-Dalsgaard, J., Hen- 
ning, R., et al. 1998, In: Provost J., Schmider 
F.X. (eds.) IAU Symp. 181: Sounding Solar and 
Stellar Interior (poster volume). OCA & UNSA, 
Nice, p. 161 

Charbonnier, P., Blanc-Feraud, L., Aubert, G., Bar- 
laud, M. 1997, IEEE Trans, on Image Processing 
6(2), 298 

Corbard, T. 1997, In: Proceedings of IX IRIS meet- 
ing, in press 

Corbard, T., Berthomieu, G., Morel, P., et al. 1997, 
A&A 324, 298 

Corbard, T., Berthomieu, G., Provost, J., Morel, P. 
1998a, A&A 330, 1149 

Corbard, T., Blanc-Feraud, L., Berthomieu, G., 
Provost, J. 1998b, A&A, to be submitted 

Duvall Jr, T.L., Dziembowski, W.A., Goode, P.R., 
et al. 1984, Nature 310, 22 

Elliot, J.R. 1997, A&A 327, 1222 

Gough, D.O., Sckii, T. 1998, In: Provost J., Schmider 
F.X. (eds.) IAU Symp. 181: Sounding Solar and 
Stellar interior (poster volume). OCA & UNSA, 
Nice, p. 93 

Gough, D.O., Thompson, M.J. 1991, The in- 
verse problem. In: Cox A.N., Livingston W.C., 
Matthews M. (eds.) Solar Interior and Atmo- 
sphere. The University of Arizona Press, Tucson, 
p. 519 



Hadamard, J. 1923, Lectures on the Cauchy Prob- 
lem in linear Partial Differential Equation, Yale 
University Press, New Haven 

Hansen, C.J., Cox, J.P., Van-Horn, H.M. 1977, ApJ 
217, 151 

Hansen, P.C. 1992, SIAM Review 34, 561 

Hansen, P.C, Mosegaard K. 1996, Numerical Linear 
Algebra with Applications 3(6), 513 

Hansen, P.C, Sekii, T., Sibahashi, H. 1992, SIAM J. 
Sci. Stat. Comput. 13, 1142 

Kirsch, A. 1996, An Introduction to the Mathemat- 
ical Theory of Inverse Problems, Springer- Verlag, 
New York 

Kosovichev, A.G. 1996, ApJ 469, L61 

Kosovichev, A.G. 1998, In: Provost J., Schmider 
F.X. (eds.) IAU Symp. 181: Sounding Solar and 
Stellar Interior (poster volume). OCA & UNSA, 
Nice, p. 97 

Morel, P., Provost, J., Berthomieu, G. 1997, A&A 
327, 349 

Schou, J., et al. 1998, ApJ, submitted 

Spiegel, E.A., Zahn, J.P. 1992, A&A 265, 106 

Teboul, S., Blanc-Feraud, L., Aubert, G., Barlaud 
M. 1998, IEEE Trans, on Image Processing 7(3) 

Thompson, M.J. et al. 1996, Science 272, 1300 

Tikhonov, A.N. 1963, Sov. Maths. Doklady 4, 1035 

Tomczyk, S., et al. 1995, Solar Phys. 159, 1 

Vogel, C.R., Oman, M.E. 1996, SIAM Journal of Sci- 
entific Computing 17(1) 



